program number2
use analyse_lattice
integer no,np
type(real_8) y(nvd)   ! nvd=3
type(normalform) normal
type(damap) map
type(taylor) x2ave
real(dp) fix(nvd),prec,xran,xran0,br(nmult),bs(nmult)
integer i,pos,k

pos=1

call make_lattice("lattice36.txt")

write(6,*) " random seed "
read(5,*) xran
write(6,*) "  width of each random multipole "
read(5,*) (br(i),i=1,nmult)
write(6,*) "  value of systematic multipole "
read(5,*) (bs(i),i=1,nmult)
XRAN0=XRAN

do i=1,lattice%n
 if(lattice%magnet(i)%type/=0) then
   do k=1,nmult
    lattice%magnet(i)%k(k)=lattice%magnet(i)%k(k)+ (2*bran(xran)-1d0)*br(k) + bs(k)
   enddo
 endif
enddo

open(unit=20,file="results_x2ave.txt")

! Compute One-turn map and normalize it including Chromaticities
fix=0.d0;


call find_orbit_with_fpp(lattice,pos,fix)

no=4;np=1; pos=1

call init(no,1,np,0)
call alloc(y); call alloc(normal); call alloc(map);


do i=1,nvd
y(i)=fix(i)+(1.d0.mono.i)
enddo

call track(y,lattice,pos)

do i=1,2
map%v(i)=y(i)
enddo
normal=map

write(20,*) " random seed "
write(20,*) xran0
write(20,*) "  width of each systematic multiple "
write(20,*) (bs(i),i=1,nmult)
write(20,*) "  width of each random multiple "
write(20,*) (br(i),i=1,nmult)

 
 write(20,*) "   "
 write(20,*) " Tunes "

call print(normal%dhdj%v(1),20,prec)

 write(20,*) "   "
 write(20,*) " Average at every Element "


! Now computes the average of 2*x**2   i.e.  <2 x**2> in terms of invariants
call alloc(x2ave)

do i=1,2
 y(i)=fix(i)+normal%a_t%v(i)
enddo
y(3)=fix(3)+(1.d0.mono.3)

prec=1.d-10

do i=1,lattice%n
 x2ave=2*y(1)**2
 call find_average_in_circle_space(x2ave)
 write(20,*) " "
 write(20,*) i,lattice%magnet(i)%twiss(-1),lattice%magnet(i)%name
 call print(x2ave,20,prec)
 
 call track(y,lattice,i,i+1)


enddo
 x2ave=2*y(1)**2
 call find_average_in_circle_space(x2ave)
  write(20,*) " "
  write(20,*) lattice%n+1,lattice%magnet(lattice%n)%twiss(0),lattice%magnet(1)%name
 call print(x2ave,20,prec)

close(20)
call kill(x2ave);




call kill(y); call kill(normal); call kill(map);

call write_lattice("lattice36_100_20.txt")


end program number2

subroutine find_average_in_circle_space(t)
use polymorphic_complextaylor
type(taylor), intent(inout):: t
type(taylorresonance) tr
external filter_ave
real(dp) filter_ave

call alloc(tr)

tr=t

call cfu(tr%cos,filter_ave,t)

call kill(tr)


end subroutine find_average_in_circle_space

  function filter_ave(j)
  use polymorphic_complextaylor
    implicit none
    real(dp) filter_ave
    integer i
    integer,dimension(:)::j

    filter_ave=1.d0
    do i=1,c_%nd
       if(j(2*i)/=j(2*i-1)) filter_ave=0.d0
    enddo

  end  function filter_ave
